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ABSTRACT 

Context. The diffusive shock acceleration mechanism has been widely accepted as the acceleration mechanism for galactic cosmic 
rays. While self-consistent hybrid simulations have shown how power-law spectra are produced, detailed information on the interplay 
of diffusive particle motion and the turbulent electromagnetic fields responsible for repeated shock crossings are still elusive. 

Aims. The framework of test-particle theory is applied to investigate the effect of diffusive shock acceleration by inspecting the 
obtained cosmic-ray energy spectra. The resulting energy spectra can be obtained this way from the particle motion and, depending 
on the prescribed turbulence model, the influence of stochastic acceleration through plasma waves can be studied. 

Methods. A numerical Monte-Carlo simulation code is extended to include collisionless shock waves. This allows one to trace the 
trajectories of test particle while they are being accelerated. In addition, the diffusion coefficients can be obtained directly from the 
particle motion, which allows for a detailed understanding of the acceleration process. 

Results. The classic result of an energy spectrum with is only reproduced for parallel shocks, while, for all other cases, the energy 
spectral index is reduced depending on the shock obliqueness. Qualitatively, this can be explained in terms of the diffusion coefficients 
in the directions that are parallel and perpendicular to the shock front. 

Key words. (ISM:) cosmic rays — plasmas — shock waves — turbulence — acceleration of particles — diffusion 


1. Introduction 

For many decades now, astrophysicists have investigated the 
origin of cosmic rays, the high-energy charged particle radi¬ 
ation that permeates space. Various experiments, both earth- 
bound and space-borne, measure a characteristic cosmic-ray en¬ 
ergy spectrum, which follows a (more-or-less) universal power 
law. The individual particles have kinetic energies above sev¬ 
eral MeV, and even energies up to more than 10^** eV have been 
detected experimentally, as shown by air-shower experiments 
(e. g., Abraham et al., 2008; Letessier-Selvon & Stanev, 2011). 
Overall, the data come from a multitude of astrophysical exper¬ 
iments that have been conducted over past decades. The under¬ 
standing of the origin of the distinct shape of the spectrum con¬ 
tinues to be an active field of study, and new experiments like 
the Square Kilometer Array (SKA, e. g., Falcke et al., 2004) are 
being assembled. 

While ineffective on its own, the original Fermi (1949) ac¬ 
celeration mechanism inspired later research, leading to the de¬ 
velopment of the so-called diffusive shock acceleration mecha¬ 
nism or of first-order Fermi acceleration (Axford et al., 1977; 
Krymsky, 1977; Blandford & Ostriker, 1978; Bell, 1978). This 
theory outlines an effective acceleration mechanism for charged 
particles at magneto-hydrodynamic shock waves. Such shock 
waves of varying sizes are ubiquitous throughout space, from 
solar winds through supernova remnants (SNR) to relativistic 
shocks in exotic cosmic objects, such as active galactic nuclei 
or pulsars (see Balogh & Treumann 2013 for an introduction). 

Thanks to the efficiency of the diffusive shock acceleration 
mechanism and the existence of plausible acceleration sites, it 
became accepted as the de-facto standard acceleration mecha¬ 
nism (cf. Abdo et al., 2010; Helder et al., 2012). There are, how¬ 


ever, still many open questions that cannot be answered entirely 
yet. These include direct proof for acceleration at the shock, the 
injection problem, and the origin of ultra-high energy cosmic 
rays has not yet been understood (e. g., Abbasi et al., 2014). In 
general, the complexity of the physical problem at hand requires 
the application of multiple approximations to find results both 
analytically and numerically. 

The work presented here will concentrate on non-relativistic 
shock waves in supernova remnants (for a review see, e. g., 
Reynolds, 2008; Blasi, 2013). Based on the numerical model 
used in the Padian code (Tautz, 2010a), a test-particle simula¬ 
tion is extended to include shock waves. This approach is then 
applied in order to investigate the effect of diffusive shock ac¬ 
celeration by inspecting the obtained cosmic-ray energy spectra. 
This relatively simple test-particle model is compared to theory, 
on one hand, and the results of Caprioli & Spitkovsky (2014), 
on the other hand, who applied a sophisticated kinetic hybrid 
model to investigate particle acceleration at non-relativistic as¬ 
trophysical shocks. Finally, the work presented here looks at the 
influence of the model of turbulence on the acceleration pro¬ 
cess. Most notably, a simple magnetostatic model is compared 
to an extended model based on Alfven waves, following the re¬ 
search of Schlickeiser (2002) on the transmission of these waves 
through a parallel shock front. 

This article is organized as follows. In Sec. 2, a brief outline 
of diffusive shock acceleration is given, followed by a summary 
of the numerical model in Sec. 3. The results are presented in 
Sec. 4. Section 5 provides a brief summary and a discussion of 
the results. 
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2. Diffusive shock acceieration 

Shock waves are a special case of a general physical phe¬ 
nomenon, where the physical properties of a (magneto-)hydro- 
dynamic environment change discontinuously across a surface 
area. They are characterized by being propagating structures; 
i. e., they have a non-zero mass flux across the discontinuity sur¬ 
face, which is commonly referred to as the shock front. In re¬ 
ality, a shock front is an irregular surface area and the discon¬ 
tinuity will occur over a certain small, but non-zero, distance. 
Any realistic description of such a system would be highly com¬ 
plex. Instead, one often approximates the shock front to be one¬ 
dimensional, planar, and to have infinitesimally small thickness. 
This simplification is justified on any scale far smaller than the 
total scale of the shock wave but much larger than the shock 
front’s thickness. 

Shock waves with very high upstream Mach numbers M\ » 
1 are called strong shock waves. In this limit, the compression 
ratio r, i. e. the ratio between the upstream density pi and the 
downstream density p 2 , can be written as (e. g., Schlickeiser, 
2002, Sec. 16.2.2.1) 

j. _ _ T + ^ ^ (1) 

Pi V 2 y-\+ i^im\ r -1 ’ 

Here, the plasma fi is the ratio of the plasma pressure to the mag¬ 
netic pressure, which will be assumed to be negligible small 
compared to the square of the shock Mach number Mj. The 
right-hand side is written in terms of the ratio of specific heats of 
the gas, y. The compression ratio reaches its minimum of r = 4 
for a monoatomic gas. An example for a strong shock wave is 
the shell of a supernova remnant (SNR), which expands with a 
velocity in the order of lO'^kms^^. Compared to the speed of 
sound, the strong shock approximation is plausible. 

Shocks with a velocity Vshock larger than the Alfven 
speed Va, which can be approximated to lie in the range of 
30 km s“^ to 50km s"' in the interstellar medium (cf. Kang, 
2013), are said to be super-Alfvenic. These fast shocks are com¬ 
mon elements of astrophysical environments. They can be found 
for example in SNRs, the relativistic bulk outflow from active 
nuclei, pulsar wind nebulae, or gamma-ray bursts. 

Averaging over particles crossing the shock at arbitrary an¬ 
gles and taking into account that, due to scattering processes 
both in the upstream and the downstream regions, the velocity 
vector will be quickly randomized, a particle’s energy gain for a 
complete round-trip is given by 

(§)=!¥■ 

Bell (1978) showed that only downstream are a small fraction 
of particles advected from the acceleration region of a non- 
relativistic strong shock. This leads to a probability for the parti¬ 
cle to remain within the acceleration region of P = 1 - Vshock/c. 
The resulting energy spectrum power-law is then given by 

dN{E) oc (3) 

Looking at SNR shells as potential sites for diffusive shock 
acceleration, Blandford & ©striker (1978) estimate that at least 
10% or even up to 50% of the energy of a SNR in the typi¬ 
cal standard Sedov solution with E - 10^' ergs can be used for 
particle acceleration. The limiting factor is rather the maximum 
size of the SNR shock front, which it reaches at the end of the 
so-called adiabatic Sedov-Taylor expansion. 


At the upper end of the cosmic-ray spectrum, particles are 
measured with such extreme energies that their gyroradii would 
exceed the dimensions of a SNR. Additionally, the time spent by 
particles in the vicinity of the shock front plays a role: The faster 
a particle becomes, the faster it will be ejected from the acceler¬ 
ation zone. Hillas (1984) researched this topic extensively in the 
context of ultra-high-energy cosmic rays. By equating the parti¬ 
cle’s gyroradius with the dimension L of an acceleration site, he 
finds a plausible upper bound to the maximum energy a particle 
can be accelerated to as 

E^^ = ZeBVsL. (4) 

This purely geometrical consideration is nowadays called the 
Hillas criterion (1984), and applies to any kind of acceleration 
site. For a single nucleon with Z = 1 and a supernova with a 
size L - 100pc and velocity Vs = lO^kms^* in a magnetic 
field of strength B - 10 '*^T, the above equation yields an up¬ 
per limit to the particle energy of around 3 x lO'^ eV. Cosmic 
rays consisting of heavier atoms with more nucleons, such as 
iron, could potentially be accelerated to even higher energies of 
around lO'^eV in the strong magnetic fields present in young 
SNR (Longair, 2011, sec. 17.5.3). Beyond that regime, either 
different acceleration mechanisms must play a role, or other cos¬ 
mic objects become the relevant acceleration sites. The Hillas 
plot (cf. Vukcevic & Schlickeiser, 2007; Shalchi et al., 2009) 
shows a straight line for the solution of the Hillas criterion for 
a proton with £max = 10^*^ eV. Of the various potential acceler¬ 
ations sites shown in the plot, only those above the line could 
possibly accomplish such an acceleration. 

Another open question is the origin of the initial particles. 
Both versions of the Fermi acceleration shown above assume 
that the induced particles already have relativistic velocities. 
Particles without a sufficiently high initial energy would not be 
accelerated due to ionization losses. This problem is called the 
injection problem. Various theories have been proposed to over¬ 
come this issue, which have been summarized in the review of 
Maikov & Diamond (2001). Of particular interest to the work 
presented here is the “thermal leakage” scenario, where the gen¬ 
eration of Alfven waves by accelerated particles is proposed. 
These waves could on one hand self-regulate the thermal par¬ 
ticle injection and on the other hand lead to efficient pitch-angle 
scattering, confining accelerating particles to the shock wave 
and thereby increasing the efhciency of the diffusive accelera¬ 
tion mechanism (Maikov, 1998; Winske & Quest, 1988). Based 
on this idea, Caprioli & Spitkovsky (2014) used a Maxwell- 
Boltzmann distributed injection profile and could reproduce an 
efficient acceleration process, leading to a power spectrum with 
a spectral index of -1.5 in conformance with the theory under 
the assumption of non-relativistic particle speeds. 

3. Numerical model 

In the Padian code, a Monte-Carlo method is applied to calcu¬ 
late the transport parameters of cosmic-ray test-particles. The 
particle equation of motion—the Newton-Lorentz equation—is 
solved numerically over a defined period of time, e. g., several 
thousand gyro periods. This procedure is repeated for a suffi¬ 
ciently large number of particles, until, the individual results are 
combined and statistically evaluated to reach an overall result. 

3.1. Magnetic field 

Padian splits the magnetic field into two components: the tur¬ 
bulent component SB and the mean magnetic field Bq, which 
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models the large-scale magnetic background field. In the work 
presented here, a uniform constant magnetic background field 
directed along the z-axis is applied. 

The turbulent component SB is calculated for every spatial 
coordinate (x, y, z) and is not discretized. A detailed description 
of the model of magnetic turbulence applied in Padian is given 
in Tautz (2010a) and Tautz & Dosch (2013). The following will 
reiterate the most important aspects of it. 

Based on the work of Shalchi & Weinhorst (2009), the power 
spectrum of the turbulence is modeled by a modified kappa dis¬ 
tribution 


G(k) = 


ihky 


( 5 ) 


where q and s define the energy range and inertial range spectral 
indices, respectively. The transition between these two ranges 
occurs at the so-called bend-over scale Iq, which will be used for 
normalization purposes as discussed below. For q - 0 and s - 
5/3 one obtains the Kolmogorov spectrum of turbulence with a 
flat energy range. An important alternative to the above spectrum 
has been given by Iroshnikov (1964) and Kraichnan (2003), who 
extended the original theory to fluids carrying a magnetic field: 
In their model, which is nowadays called Iroshnikov-Kraichnan 
turbulence, they assume locality of interactions, weak interac¬ 
tions, and isotropy. For turbulence in a plasma influenced by a 
uniform magnetic field they then arrive at an energy spectrum of 
the form 

E(k) cc (eVAy^^k-^'\ (6) 


where Va - BqI ^jAnp denotes the Alfven velocity. Numerous 
experiments successfully reproduced either form of the en¬ 
ergy spectrum when investigating turbulent systems (Bruno & 
Carbone, 2013). But both of these theories, and many others, as¬ 
sume an isotropic turbulence which is contradicted by many ob¬ 
servational results, in particular in the solar wind. More recently, 
Goldreich & Sridhar (1997) presented another model that takes 
the anisotropy in the turbulence of astrophysical magnetohydro¬ 
dynamic (MHD) plasmas into account. However, no model has 
been found which can explain all experimental data nor predict 
the results of numerical simulations. A recent review on this mat¬ 
ter is given by Schekochihin & Cowley (2007). 

Padian furthermore follows the ideas of Giacalone & Jokipii 
(1999) and generates the turbulent magnetic field by superposing 
multiple plane waves with random phase angles and directions 
of propagation. According to Batchelor (1953), one can model 
isotropic, spatially homogeneous turbulence this way, when us¬ 
ing a large number of waves modes. For an individual mode, 
the magnetic field can then be calculated with (Tautz & Dosch, 
2013) 

SBn = ixAikn) cos (kn( + l/l„) , (7) 

where A(kn) oc sjG(kn)Akn is the amplitude function, which is 
defined by the energy spectrum (cf. Eq. (5)). The unit vector 
is perpendicular to the direction of which in turn is obtained 
by applying a random rotation matrix to the spatial coordinate 
{x,y,z) for each wave mode A:„. Additionally, each wave mode 
has a random phase given by i/i„. 

The real part of a summation over a number of plane wave 
modes N then yields the total turbulent magnetic field as 


3.1.1. Alfvenic turbulence 

The turbulence model introduced above can be extended to sup¬ 
port propagating Alfven waves (e.g., Michalek & Ostrowski, 
1996; Tautz, 2010a). This is achieved by modifying the turbu¬ 
lent magnetic field components of individual plane wave modes 
SBn in Eq. (8) with an oscillation frequency, yielding 

SBn = IxAikn) cos [k„( - a>(kn)t + l//n] ■ (9) 

Here, w(k„) denotes the dispersion relation of the MHD wave, 
which, for Alfven waves, is given by w(k„) = ±Va^II- It has to 
be noted that this model is only valid in absence of any shock 
effects, i. e. in the upstream region. Sec. 3.2.1 below discusses 
the downstream region. 

Optionally, Padian can also include the effect of the electric 
field components of the Alfven waves. These are perpendicular 
to the magnetic field components and obtained from applying 
the Earaday induction equation 

B^-kxE (10) 

ti) 

by taking the polarization properties of Alfven waves into ac¬ 
count. 

3.2. Shock wave 

In order to investigate the effects of diffusive shock acceleration, 
Padian was extended to model a one-dimensional MHD shock 
front: The ambient plasma will flow through the shock rest frame 
with a defined velocity toward the shock front. At the position of 
the shock front, a discontinuity in the magnetic field and ambient 
gas velocity is added according to the Rankine-Hugoniot jump 
conditions (see, e.g., Schlickeiser 2002, Chpt. 16 and Longair 
2011, Chpt. 11). Note that the back reaction of the cosmic-ray 
particles on the shock wave will be neglected (cf. Riquelme & 
Spitkovsky, 2010). 

Padian internally works with dimensionless parameters, 
which simplifies the support for relativistic particles and im¬ 
proves the generality of the results. The time t is normalized with 
the relativistic gyrofrequency to t = = qBotHymc) with y 

the relativistic Lorentz factor. In place of the momentum, Padian 
instead uses the particle’s rigidity R normalized to the turbulence 
bend-over scale Iq, as given by R = RlI^o - vfiQieito- 

Using these parameters, the equation of motion is solved in 
the gas rest frame instead of in the observer frame. According to 
Earaday’s law, the flow of the ambient plasma will also induce 
an electric field 

£'ind. =--/?flowXB, (11) 

C 

where Rgow represents the flow rigidity, and B the total magnetic 
field. In the shock frame of reference, which is used in the simu¬ 
lations presented in Sec. 4, this will influence the particle rigid¬ 
ity even without them crossing the shock front. The effect of the 
flow of the ambient plasma gas on the test particles is taken into 
account by adapting the Newton-Lorentz equation to include the 
inductive electric field from Eq. (11), yielding 

-i/? = (R - Rflow) X B. (12) 

dr 


N 

SB ^%Yj5Bn. 

n=l 


The relative rigidity Rfiow follows from the Lorentz transforma- 
(8) tion into the gas rest frame. It is always parallel to the shock nor¬ 
mal, and scales analogously to the gas velocity across the shock 
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front. For a given value of the shocks rigidity /?s and compres¬ 
sion ratio r, one can calculate the relative rigidity upstream and 
downstream of the shock front using 

^rel,l = Rs, (13a) 

/?rel,2 = -RS (13b) 

r 

Upstream, the magnetic field B\ is calculated directly according 
to Eq. (8). Downstream, the tangential magnetic field compo¬ 
nent B 2 ,t is scaled with the compression ratio as defined by the 
Rankine-Hugoniot equations, while the normal component is 
continuous, leading to 

B 2 ,n — Bn (14a) 

B2,t = rB,. (14b) 

This simple model allows the investigation of parallel, oblique, 
and perpendicular fast, strong shock waves. 


treatment holds only for unidirectional circularly polarized up¬ 
stream wave fields, since variations in the magnitude of the per¬ 
pendicular magnetic field component may induce motion of the 
shock front itself from its equilibrium position (e. g., Achterberg 
& Blandford, 1986). The situation gets even more complicated, 
if the upstream wave vectors are not aligned with the shock nor¬ 
mal, since the assumption of a planar shock is then probably 
violated, also. This is why we do not expect our model to be 
generalized to oblique shocks very easily.” 

3.2.2. Spatial diffusion 

In presence of a shock, the typical way of computing the spatial 
diffusion coefficients, /fx.y.z, must be altered to take the ambient 
gas flow velocity into account. Originally, these coefficients are 
obtained from the particles’ mean square displacements via 

^/.orig ~ T 2 ^(“^/(tmax) “ 2C/(0)) ^ , 1 E {x, y, z}. (18) 

4-»max 


3.2.1. Interaction with Alfven Waves 


Alfven waves interact with parallel shock fronts, which process 
was described in detail by Campeanu & Schlickeiser (1992); 
Schlickeiser et al. (1993); Vainio & Schlickeiser (1998, 1999), 
and Vainio et al. (2003). A review of this combined work is given 
in Schlickeiser (2002, Chpt. 16.3). 

The important result for the work presented here are the 
transmission coefficients at constant wavenumber k. They are 
found to be (Schlickeiser, 2002, eq. 16.3.26f) 


^^r+ 1 / M + We.i 
2^J? \ M + Vr//c,i/ 
V7- 1 / M + //e.i 
2V? \ M - 2/rHc^il 


(15a) 

(15b) 


In the shock rest frame, the gas flow upstream and downstream 
of the shock front would also be included in the equation above. 
This is undesired, as one is rather interested in the diffusion 
perpendicular or parallel to the magnetic field lines, which, in 
the flux-freezing approximation, flow together with the ambient 
plasma. To account for that, time-dependent “running” diffusion 
coefficients are used (Tautz & Shalchi, 2010), which are flow- 
corrected and read 

Ki(t) = Yf ~ -^'(O) “ 4xflow,i(0)^) ■ (19) 

Here, the correction factor zlji:flow(0 is introduced, which corre¬ 
sponds to the displacement of the initial particle location Jto due 
to the ambient gas flow. It is defined as 


Here, r is the shock compression ratio and s is the inertial range 
spectral index with values 1 < s < 2. For the work presented 
here, it is set to 5/3 (cf. Sec. 3.1). Finally, the equation above in¬ 
cludes the normalized cross helicity state of the waves, which 
is given by 


(6Bff - (^Rb) 
((5Bf)2 + (dZ?b)2' 


4lA^flow(0 “ f RnowC^oC^ )) tit , (20) 

Jo 

where /?flow(Jto(0) is either the upstream or downstream flow ve¬ 
locity, depending on the relative position of Ji:o(f) to the shock 
front. In the Radian simulation, this factor can only be approxi¬ 
mated by a discrete summation over a number N of time steps 


In Radian, only forward moving waves are injected upstream so 
that ffc,i = +1. There, the turbulent magnetic field is calculated 
using Eq. (9). Downstream of the shock front, the formula is 
adapted to take both the forward and backward traveling Alfven 
waves with the corresponding transmission coefficient into ac¬ 
count, leading to 

(5^2 = ^ ^ [Tk sin (k„( - uik„)t - 'P„) 

n=l 

+ Rk sin {k„^ + (Aj(k„)t - Wn)] ■ (17) 

This model for the downstream turbulent magnetic field already 
leads to a significant increase in the magnetic field strength, 
compared to the upstream values. As such, contrary to the mag¬ 
netostatic model of turbulence of Eq. (8), no artificial compres¬ 
sion following Eq. (14b) must be applied. 

Finally, it is important to state that this model is currently 
limited to parallel shock configurations. Vainio & Schlickeiser 
(1999) explicitly mention in their work that it might not be eas¬ 
ily extensible to oblique shock configurations: “Note that our 


N 

^•^flow(0 “ ^ ] ^flow(-t:Q(u zlt))zlt, (21) 

H=0 

where NAt — t. It is expected that the impact on the final values 
of the diffusion coefficient is negligible because the flow velocity 
is small compared to the particle velocity. 

3.3. Injection 

The work presented here investigates two methods of parti¬ 
cle injection. In the simplest case, all particles of the simula¬ 
tion are injected with a constant absolute rigidity value. The 
second injection profile is inspired by the work of Caprioli & 
Spitkovsky (2014) and follows the thermal leakage model out¬ 
lined by Maikov & Diamond (2001). There, the particles are ini¬ 
tially assigned an absolute non-relativistic rigidity value follow¬ 
ing an offsetted Maxwell-Boltzmann distribution of the form 

R - RoSs + fc amma (/?;3/2,/3). 


( 22 ) 
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Rigidity R 


yield 


d / 

AnAE— ^N(E,AE)v, 
aE 


(26) 


where v is the particle’s velocity. The above equation can be 
transformed to use the normalized particle rigidity R applied by 
Padian, instead of the energy. Starting with the energy 


E - jmc^. 


(27) 


one can substitute the Lorentz factor y = -^1 -i- (p/mc)^ yielding 


E = + (pcf. (28) 

Replacing the momentum with p - RQmlo and simplifying the 
expression leads to 


Fig. 1. Rigidity injection profile following a Maxwell distribution ac¬ 
cording to (22) with p = (Rmax - Rmm)/2 and Roffs = Rmm, with 
f^min = 10“'' and = 10“^ *. The red line corresponds to lO"* particles 
used in a simulation, obtained from the original Maxwell distribution 
shown by the black dotted line with a cutoff value of Rcutoff = 10“^. 


E = Qmkc H- (29) 

with Rc - cIQIq. Under the assumption of R » R^, which is 
valid as long as v/c » 1 / V2, one arrives at E cc R. The equation 
for the power spectrum can thus be written as 


Here, a library function is used for the gamma distribution. Its 
probability density function in terms of the shape factor a and 
rate factor fi is given by 

/Gamma(Ma,j6) = , (23) 

Eia) 

where E denotes the gamma function. For a shape factor a - 
3/2 this yields r{a) - •0r/2. Substituting/I = Ijk^T, one then 
obtains the Maxwell-Boltzmann distribution function of energy 


d; N{R,AR) 
d£ RAR 


(30) 


In the non-relativistic limit R « /?c on the other hand, Fkin 
leading to 


dj N(R^,AR^) 

- OC ----- 

dFkin R^AR^ 


(31) 


In the spectra plots below, logarithmic spacing is chosen by sam¬ 
pling the above equation with a constant ratio ARfR. 


rp / 1 

/mw(£) = 2 y- J (24) 

The rate parameter /3 is configurable in the simulations and de¬ 
fines the width of the distribution. The offset value Rofts intro¬ 
duced in Eq. (22) moves the whole distribution into a user de¬ 
fined area. 

By default, this distribution would also produce very low 
rigidity values, which are uninteresting for the shock accelera¬ 
tion process, as their energy will not change considerably. Thus, 
to improve the runtime performance of the simulation, the dis¬ 
tribution function is additionally modified to include a cutoff at 
^cutoff- Rigidity values below this value are discarded and the 
distribution function queried again. Combined, this scheme pro¬ 
duces an injection profile which is shown in Fig. 1. 

In both injection profiles, the initial rigidity direction of the 
individual particles is randomized. Additionally, a distinct tur¬ 
bulent magnetic field is created for each particle by varying the 
seeds of the random number generators used in the calculation 
of Eq. (8). 

3.4. Acceleration spectrum 

An energy spectrum comparable to the measured cosmic-ray en¬ 
ergy spectrum can be obtained from the numerical results as fol¬ 
lows. Eollowing the ideas of Giacalone & Jokipii (1996), the 
number of particles with a momentum in the range [p, p+Ap^is 
given by 

N{p,Ap) = 4np^fip)Ap, (25) 

where f(p) represents the distribution function for the momen¬ 
tum p. In terms of the flux dj/dE this can then be rewritten to 


4. Simulation results 

The following chapter presents various results obtained from 
running Padian under different parameter configurations. To re¬ 
duce the parameter space, the following values where set in all 
simulations: 

- The compression ratio r is set to 4. 

- The magnetic mean field Bq is set to be uniform and points 
along the z-axis. 

- Two types of turbulent magnetic field 6B generation are 
studied. On one hand, the magnetostatic case according to 
Eq. (8), or via Alfven waves and Eq. (9). In both cases, 
Am = 128 wave modes are superposed in slab geometry. 
A Kolmogorov energy spectrum is chosen for the magnetic 
turbulence, i. e. i = 5/3, with an additional energy range pro¬ 
portional to kP, here with q -Q. The minimum wavenumber 
exponent in units of the bend-over scale is set to -5 and the 
maximum one to 3. Finally, the turbulent magnetic field is 
normalized to fulfill the condition SB IBq - 1, corresponding 
to a model of strong turbulence. 

- For every particle injected into the simulation, a separate tur¬ 
bulent field is created with different random seed values. This 
ensures the independence of individual test particles when 
analyzing all as a statistical ensemble. 

For all simulations, the particle energy spectrum at differ¬ 
ent times is obtained. This allows one to estimate: (i) the in¬ 
fluence of the initial particle energy spectrum; (ii) the method 
of magnetic held turbulence generation; and (iii) the direc¬ 
tion of the shock has on the final spectrum. Within the result¬ 
ing particle energy spectra, a range is selected and its double- 
logarithmic data fitted linearly to And the spectral index. Note 
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Time T 

Fig. 2. Parallel magnetostatic shock with particles injected at a fixed 
rigidity of R(0) = 10“^. The time evolution of the power spectrum, i. e., 
the flux dj/dR over total rigidity gain R/R(0), over the runtime of a 
Padian simulation of a strong parallel shock with 10“* test particles is 
shown. The feature on the left side indicates that some particles lose 
energy. To the right, a smooth constant power law establishes itself over 
time. 



Relative Rigidity RIRo 


Fig. 3. Same as Fig. 2 but only for the last time step. A linear fit of the 
double logarithmic data of the final power spectrum for the last time 
step at T = 10® is shown. For the data range spanned by the blue fit line, 
a spectral index x » -2.01 ± 0.03 is obtained. 


that the errors given for the spectral index come from this analy¬ 
sis alone. While all results are qualitatively reproducible, differ¬ 
ently seeded Padian simulations with otherwise equal configura¬ 
tion parameters would yield slightly different numerical results 
in the order of the fit errors. 

4.1. Benchmark results 

Initially, different configurations of a system with magnetostatic 
turbulence are investigated. The turbulent magnetic field com¬ 
ponent is generated using Eq. (8) and can easily be adapted to 
oblique shocks. 

First, the effects of shocks on an ensemble of 10^ particles 
injected at the origin with an equal rigidity value of R = 10 “^ 
is investigated. All particles are assigned random initial rigidity 
directions. The shock rigidity is set to Rs = f?/ V3 5.8 x 10“^, 
following the previous work by Giacalone & Jokipii (1999). The 
simulation runs in the shock rest frame, where the shock front 
lies in a plane through the origin, such that all particles are in¬ 
jected directly at the shock. This system is then simulated until 



Rigidity R^ 


Oe 00 2e 04 4e 04 6e - 04 8e 04 le 05 
Time r 

Fig. 4. Parallel magnetostatic shock with a Maxwellian injection profile 
as indicated by the black dotted line. This figure shows how the ini¬ 
tial power spectrum is modified by the interaction with the shock front. 
Particles get accelerated and, over time, form a smooth constant power 
law, close to parallel to the abscissa. 



Fig. 5. Same as Fig. 4 but only for the last time step. A linear fit of the 
double logarithmic data of the spectrum for the last time step at r = 10® 
yields a spectral index of -1.49 ± 0.02. 

the normalized time t = Q^eit - 10® is reached with the rela¬ 
tivistic gyrofrequency. At 100 equidistant time points, the power 
spectrum of all particles is computed according to Eq. (30). 

For a parallel shock, where the shock front lies in the x-y- 
plane, the results are shown in Fig. 2. One can observe multi¬ 
ple effects; First, note how a fraction of the particles lose en¬ 
ergy. This effect is a result of the relative motion of the ambient 
gas on the rigidity of the injected particles and unrelated to the 
shock acceleration. Second, all other particles are efficiently ac¬ 
celerated by the shock. Already after the first snapshot, a power 
law is exhibited for parts of the energy spectrum. Over time, 
this range grows as particles continue to be accelerated. A lin¬ 
ear fit of the double-logarithmic power spectrum data in this 
range, as is shown in Fig. 3, yields a spectral index of about 
X -2.01 + 0.03. This value is compatible with the theoreti¬ 
cally predicated value of -2. 

Inspection of the individual particles’ relative rigidity boost 
over the shock distance shows some interesting insights. First, 
some particles are still close to the shock front and thus within 
the acceleration zone, which explains the continuous rise of the 
total energy of the particle ensemble. Second, it appears that the 
bulk of particles gets advected downstream of the shock. These 
particles have a distance close to -Rsr/r, indicating that, un- 
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Oe-^OO 2e 04 4e-f04 6e 04 8e 04 le 05 


Time t 

Fig. 6. Perpendicular magnetostatic shock with a Maxwellian injection 
profile as indicated by the black dotted line. Again, the injection profile 
is modified by the shock acceleration, but the slope is not parallel to the 
abscissa. 



Fig. 7. Same as Fig. 6 but only for the last time step. The double log¬ 
arithmic data of the spectrum for the last time step at t = 10^ is scat¬ 
tered more strongly around the linear fit and yields a spectral index of 
X* -1.98 ±0.04. 


der the flux-freezing assumption, they follow the magnetic held 
lines which flow with the downstream ambient plasma velocity 
of Rs/r. 

The total kinetic energy of the system increases by a fac¬ 
tor of close to 7 over the simulation time. The total aggregated 
rigidity of all particles still rises toward the end of the simula¬ 
tion. This is an indication that the simulation did not yet reach 
a steady state. By increasing the simulation time further, even 
higher maximum particle energies are to be expected. However, 
no qualitative changes to the power law are expected, just the 
high-energy tail will become smoother. 

4.2. Maxwellian injection profile 

Inspired by the work of Caprioli & Spitkovsky (2014), the effect 
of the shock acceleration on a particle ensemble with an initially 
Maxwell-distributed rigidity spectrum is also investigated. The 
simulation was adapted such that each of the 10^ particles is 
injected with an absolute rigidity value given by Eq. (22). For 
the simulations below, the distribution parameters were set to 
^ — (f^max “ f^min)/2 and feoffs — with /?min — 10 and 
f^max = 10“^ '. The rigidity cutoff was set to /^cutoff = 10“^. The 
initial direction of each particle is, as before, still random. The 


shock rigidity is set to Rs - 10“^ and only a limited time period 
up to T = 10^ is simulated. Increasing this value improves the 
smoothness of the power spectrum in the high-rigidity regime, 
at the cost of significantly longer simulation run times, but does 
not influence the results for the power spectra qualitatively. 

For comparison purposes with the original work of Caprioli 
& Spitkovsky (2014), the power spectra are now obtained by 
plotting the flux to the power of 1.5 over the square of the 
rigidity. This relation corresponds to a particles’ kinetic en¬ 
ergy dependence on the rigidity in the non-relativistic limit (cf. 
Sec. 3.4). This is necessary because the theory of the Maxwell- 
Boltzman distribution is only applicable to non-relativistic ener¬ 
gies. Contrary to the previous injection model with relativistic 
constant rigidity, one now expects spectral indices of 1.5 from 
theory (Caprioli & Spitkovsky, 2014). 

Furthermore, one cannot directly calculate diffusion coeffi¬ 
cients for this model of injection. The formula in Eq. (19) is only 
meaningful for particles of similar starting energy. Here, one 
would need to group the injected particles by energy and calcu¬ 
late individual coefficients, which requires a significant increase 
in the number of injected particles. Due to the higher computa¬ 
tional effort this would require, the spatial diffusion coefficients 
are omitted below. 

4.2.1. Parallel shock 

For a parallel shock configuration a smooth power spectrum with 
a spectral index of about -1.5 quickly establishes itself, as can 
be seen in Fig. 4. The high-rigidity tail of the power spectrum 
is still in flux and expected to smooth out for longer simulation 
times. 

It is notable how the injection profile is influenced also in the 
low-rigidity region below the shock rigidity. This effect might be 
due to the cutoff applied to the injection profile. However, while 
removing this may influence the shape of the power spectra, no 
qualitative change to the spectral index in the constant power law 
region is expected. Rather, one can expect that the direct injec¬ 
tion at the shock front leads to this result. If significantly more 
particles are injected, in a larger space around the shock front, 
the low-energy tail of the power spectrum should stay rather un¬ 
affected. But, since this requires significantly higher computa¬ 
tional effort, this assumption was not verified in the work pre¬ 
sented here. 

While Caprioli & Spitkovsky (2014, fig. 1) do not list a value 
for the spectral index they obtain for a parallel shock, the power 
spectrum also exhibits a region close to parallel to the abscissa, 
in agreement with the results presented here. Compared to the 
previous results obtained from the simulation of a constant in¬ 
jection profile discussed in Sec. 4.1, this result only differs due 
to the non-relativistic calculation of the power spectrum. 

4.2.2. Perpendicular shock 

For a perpendicular shock, the simulation was set up as before 
in Sec. 4.2.1, but the shock front now lies in the v-z-plane. The 
background magnetic held Bq still points in the z-direction, and 
is thus transverse to the shock normal. Downstream of the shock, 
the transverse components of the combined magnetic held B are 
amplified by the shock compression ratio. In this case, the results 
differ strongly from those of a parallel shock. Most notably, the 
particles are only accelerated during a short time period, as seen 
in Fig. 6. This is in agreement with previous studies (Giacalone, 
2005; Caprioli & Spitkovsky, 2014, and references therein), who 
























































8 


Wolff & Tautz: Cosmic-ray acceleration at astrophysical shocks 




time T 


Fig. 8. Flow-corrected spatial diffusion coefficients k according to 
Eq. (19) for the z-direction parallel to the magnetic background field 
Bo, and the perpendicular components x and y. For a parallel shock as 
shown in the upper panel, the diffusion in z-direction, parallel to the 
shock normal, is the highest. The perpendicular diffusion coefficients in 
X and y direction are approximately equivalent, as expected. Throughout 
the shock simulation, the values of all diffusion coefficients continue to 
rise. A perpendicular shock configuration on the other hand produces 
significantly different results, as shown in the lower panel. Here, all 
three components differ significantly from another, with the perpendic¬ 
ular y-component, parallel to the shock normal, having the lowest dif¬ 
fusion coefficient. The parallel z-component, here perpendicular to the 
shock normal, is again the highest. The r-component lies in between 
these two. Furthermore, the diffusion coefficients and /Cy show a peak 
at around r =» 250. The parallel diffusion coefficient on the other hand 
continues to rise after this time, but significantly slower than before. 


showed that perpendicular shocks accelerate particles at a higher 
rate than parallel shocks. In particular, Giacalone (2005) found 
that the flux at the highest energies is dominated by particles 
accelerated at a perpendicular shock. 

The particles diffusion coefficient in these directions is sig¬ 
nificantly larger than in the normal, y-direction, as shown in 
Fig. 8 for a constant injection profile. Once the particles cross 
the shock front, they are advected from the acceleration zone 
with the ambient downstream plasma flow. Because of that, a 
much shorter simulation time of only t = 3 x 10^ is sufficient 
to capture the dynamics of the acceleration process, indicated 
also by the plateau reached in the total rigidity. Within this time, 
the diffusion coefficients leave the initial ballistic regime (see 
Fig. 8). After reaching a peak at around t = 250, the transverse 
diffusion coefficients /c^ y decrease over time. Due to (B^ y) = 0 
one expects these values to approach zero for even higher times. 
The diffusion along z, i. e. parallel to the magnetic mean held but 


Fig. 9. Oblique magnetostatic shock with a Maxwellian injection pro¬ 
file. The spectral indices obtained for Padian simulations of strong 
shocks with varying obliqueness and a Maxwell rigidity injection pro¬ 
file. 



shock angle [°] 

Fig. 10. Spatial diffusion coefficients for various shock obliqueness an¬ 
gles, showing a variation by multiple orders of magnitude for larger 
shock angles. They seem to be correlated to the spectral indices (cf. 
Fig. 9, with higher diffusion coefficients also yielding higher spectral 
indices. Note that for the calculation of diffusion coefficients, a constant 
injection profile has been assumed. 

transverse to the shock normal, also decreases its growth at this 
time, but still increases slightly. With (B^) = Bq, one can expect 
Kz to approach a steady state for higher simulation times. The 
strong initial increase of the diffusion coefficients, especially in 
the y-direction, correlates to the acceleration process; with in¬ 
creased energy, the Larmor radii of the particles increases and 
thus the diffusion coefficients increase. 

The resulting power spectrum for this shock configuration 
is shown in Fig. 7. Fitting the more strongly scattered numer¬ 
ical data yields a spectral index with a value of -1.98 + 0.04 
corresponding to a constant power law. Note that this result is 
equivalent to that of a perpendicular shock with a constant injec¬ 
tion profile. One can still observe efficient acceleration during 
the short acceleration period, with individual particle rigidities 
being boosted by about two orders of magnitude. Overall, the 
combined system increases its kinetic energy by a factor of close 
to 2.5. 

4.2.3. Oblique shocks 

Between the two extreme cases of a parallel or perpendicular 
shock, the so called oblique shocks can be set up (Sironi & 
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Time T 

Fig. 11. Parallel Alfvenic shock with particles injected at a fixed rigidity 
of i?(0) = 10“^. The time evolution of the power spectrum shows shows 
how efficient the particles get accelerated, and that a constant power law 
steady state is approached. 



Fig. 12. Same as Fig. 11 hut only for the last time step. A (partially) 
linear fit of the double logarithmic data at the last time step is shown, 
which yields a spectral index of x x -1.95 ± 0.03. 


Spitkovsky, 2009). For simulations using a varying level of the 
shock obliqueness, the resulting spectral indices are plotted in 
Fig. 9, showing a clear angular dependency. Again, the behavior 
is qualitatively equivalent to the results obtained from the sim¬ 
ulations on a particle system with a constant rigidity injection 
profile. 

The quasi-parallel shocks produce a power spectrum with 
a spectral index slightly below -1.5. For quasi-perpendicular 
shocks at 45° and higher, the acceleration process becomes less 
efficient. The spectral index decreases significantly, which is cor¬ 
related to the less efficient acceleration at such shocks. For the 
perpendicular shocks then, a spectral index of about -2.0 is mea¬ 
sured. 

Our results are, to some extend, corroborated by the findings 
of Caprioli & Spitkovsky (2014, their Fig. 12), although their 
explanation relies on the central fact that the back-reaction of 
the particles on the shock structure is included. Accordingly, the 
structure of the self-generated turbulent magnetic fields in the 
upstream regime depends strongly on the shock angle. In con¬ 
trast, Giacalone (2005) found that, in the high-energy end of the 
particle spectrum, the spectral index tends to be compatible with 
the predictions of diffusion shock acceleration model, indepen¬ 
dently of the shock normal angle. However, it should be noted 


that they obtained the high-energy tail of the spectrum with the 
help of a special steady-state method. At intermediate energies, 
the spectral index becomes steeper with an increasing shock an¬ 
gle, similar to our findings. 

A possible explanation for our results can be found by mod¬ 
ifying the core assumptions of diffusive shock acceleration. In 
particular, the escape probability might be related to the diffusive 
behavior of the particles and, thus, might be energy-dependent. 
To illustrate this further, note that, to lowest order, particles can 
be assumed to follow the magnetic field lines, which, on aver¬ 
age, are given by the magnetic mean field Bq. The diffusion co¬ 
efficient parallel to the direction of Bq, k\\ is significantly higher 
than that in the perpendicular direction, /Cj^, as shown in Fig. 8. 
In addition, the diffusion coefficients also have a different en¬ 
ergy dependence. In quasi-perpendicular shocks, the magnetic 
field lines are advected with the downstream plasma flow, and 
the particles will then follow this flow, with Ky being the smallest 
of the three diffusion coefficient components. Additionally, the 
increased magnetic field strength downstream reduces the parti¬ 
cles’ Larmor radius. Thus the acceleration process only occurs 
for a short time period, and just a few particles reach high ener¬ 
gies, leading to the steep power law. The results for the diffusion 
coefficients for varying shock obliqueness, as shown in Fig. 10, 
also seem to validate this hypothesis. There, a correlation of the 
spectral index and the diffusion coefficients is visible: the higher 
diffusion coefficients in the direction parallel to the shock front 
enable more shock front crossings, leading to higher particle ac¬ 
celeration and thus a flatter power spectrum. For further discus¬ 
sion, the reader is referred to the investigations of shock drift 
acceleration (e. g.. Ball & Melrose, 2001; Park et ah, 2013, and 
references therein), which underline the complex interactions of 
particles with quasi-perpendicular shocks. 

4.3. Alfvenic turbulent magnetic field 

The previous results were all obtained from simulations with 
magnetostatic turbulent magnetic fields in slab geometry. In this 
section, the turbulence is generated by the Alfven-wave model 
outlined in Sec. 3.1.1 and 3.2.1, and its influence on the accel¬ 
eration process investigated. This model can only be applied to 
parallel shock configurations. In the first simulation, only the 
magnetic component of Alfven plasma waves is enabled, fol¬ 
lowed by a second simulation which includes the electric field 
component to investigate the effect of stochastic acceleration. 
All other parameters for the simulations presented in the follow¬ 
ing are copied from Sec. 4.1. Most notably, particles are injected 
with a constant rigidity of 10“^, and the shock rigidity is set to 

10-2/ V3. 

Running Padian with such a configuration yields power spec¬ 
tra as shown in Fig. 11, which again indicate a clear accel¬ 
eration process occurring at the shock front. Compared to the 
results for magnetostatic turbulence, a slightly lower value of 
X ^ -1.95+ 0.03 is measured for the spectral index of the power 
spectrum at the last time step (Fig. 12). Otherwise, no fundamen¬ 
tal differences are apparent. Likewise, the growth of the total 
kinetic energy of the particle system, as shown in Fig. 13, re¬ 
sembles the previous data for magnetostatic turbulence that was 
presented in subsection 4.2.1. 

4.3.1. Influence of the electric field 

When the electric field component of the Alfven plasma waves 
are included in the simulation, the results are significantly mod- 
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Fig. 13. Parallel Alfvenic shock with particles injected at a fixed rigidity of f?(0) = 10“^. Shown are the cases without (left panels) and with (right 
panels) turbulent electric fields. As shown in the upper panels, a strong shock efficiently accelerates the particle ensemble. The acceleration is 
strongest at the beginning, when all particles are close to the shock front. Due to parallel diffusion, many particles will not encounter the shock 
again after the initial boost, whereas others continue to get accelerated. This is apparent from the lower panels, which also show that some particles 
initially lose energy or never cross the shock to get accelerated. 


ified. Fig. 14 shows the power spectra obtained from such a 
configuration, and Fig. 13 depicts the acceleration process in 
terms of the kinetic energy boost of the total particle system as 
well as for a selection of individual particles. Compared to the 
simulation without electric fields, one sees a larger fraction of 
particles losing energy, indicated by the relatively large flux of 
particles with relative rigidities below 1 in the power spectrum. 
Furthermore, no constant power law is visible, but rather a sec¬ 
ond peak at around ~ 15 grows over time. A fit of the final 
power spectrum (Fig. 15) to find a spectral index for comparison 
reasons yields a value ofxa; -1.9 + 0.1. The total rigidity of the 
particle system (cf. Fig. 13) keeps growing after the characteris¬ 
tic strong initial boost due to shock acceleration. 

Interestingly, the total rigidity boost is of the order of one 
magnitude, whereas individual particle boost lie in the range 
of up to three orders of magnitude. The total particle ensem¬ 
ble gains energy by a factor of about 7, whereas individual 
particles encounter rigidity boosts of over more than two or¬ 
ders of magnitude. When the electric field component of the 
Alfven model of turbulence is enabled, the total rigidity boost 
for short times resembles the pure shock acceleration. Over 
time though, the influence of the stochastic acceleration due to 
the electric field becomes apparent, with the total rigidity ap¬ 
proaching a constant linear growth instead of converging to¬ 
ward a maximum value when the particles are advected away 
from the shock front (cf. Tautz et ah, 2013). In the absence of a 
shock, the stochastic acceleration due to turbulent electric fields 


is well known and has been investigated theoretically in terms of 
momentum diffusion (e.g., Skilling, 1975; Schlickeiser, 1989, 
1994) as well as numerically (e. g., Michalek et ah, 1999; Tautz, 
2010b). Time-dependent turbulence in the form of magnetohy¬ 
drodynamic plasma waves such as Alfven waves can be inter¬ 
preted as an acceleration mechanism comparable to a first-order 
Fermi process (Schlickeiser, 2009). 

In combination, the rigidity boosts of individual particles 
show the energy gain in steps, characteristic for shock acceler¬ 
ation, but also constant growths due to stochastic acceleration. 
Compared to the previous results which neglected the electric 
field, the effect of stochastic acceleration are thus clearly appar¬ 
ent. 

Note that Padian’s runtime performance significantly de¬ 
grades for high-energy particles because the numeric integration 
algorithm applies a constant error tolerance, leading to an in¬ 
crease in integration steps with growing energy. Due to that, only 
a reduced simulation time up to t = 2 x 10^ was simulated here, 
as the combined effect of both, stochastic and diffusive particle 
acceleration create too many high-energy particles. 

5. Discussion and conciusion 

In this article, the process of diffusive shock acceleration was 
investigated using an ensemble of independent test particles. In 
comparison to particle-in-cell simulations, the main advantage is 
that, aside from the model being easily adjustable, it is possible 
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Fig. 14. Parallel Alfvenic shock with particles injected at a fixed rigidity 
of R(0) = 10“^. In contrast to Fig. 11 and 12, the turbulent electric 
field components for an Alfven rigidity of = 10“* is also included. 
The figure shows how initially the acceleration process yields similar 
results to the simulation without electric fields. However, over time, 
the stochastic acceleration process influences the form of the spectra 
heavily, with a flux peak forming at around R/Ro ~ 15. 



Fig. 15. Same as Fig. 14 but only for the last time step. For compari¬ 
son reasons, a region of the final power spectrum is fitted linearly. This 
yields a qualitative result of sj -1.9 ± 0.1. 


to correlate the acceleration process to the spatial difFusivity of 
the injected particles. 

When the energy of the total particle ensemble is investi¬ 
gated, one finds power spectra being formed over time in all sim¬ 
ulations. Independent of which model turbulence was applied, 
i. e. for both the magnetostatic and the Alfven-wave model of 
turbulence, a spectral index of close to -2 was found for parallel 
shock configurations. (For non-relativistic particles, the spectral 
index is accordingly reduced to -1.5.) This value equals the one 
predicted by the theory of first-order Fermi acceleration. In both 
simulations, particles remain within the acceleration region over 
long times and the total energy of the particle ensemble keeps 
growing. However, the power spectrum establishes itself quickly 
and only the high-energy tail changes toward the end of the sim¬ 
ulation. The diffusivity of the particles over time indicates that 
the ballistic regime is not left within the maximum simulation 
time. This is due to the continued acceleration of some particles. 
To exclude this effect one either needs to apply even higher simu¬ 
lation times, or manually exclude particles from the acceleration 
region. 


The simulation with a full model of Alfven waves, including 
their electric field components, illuminates the interplay of dif¬ 
fusive shock acceleration and stochastic acceleration. This has 
some minor influences on the resulting energy spectrum, lead¬ 
ing to higher scattering of the data around the numerical fit 
corresponding to a spectral index of, in this case, -1.9 + 0.1. 
Qualitatively however, no large differences are otherwise appar¬ 
ent. 

The constructed magnetic fields of the two applied models of 
turbulence differ in the magnitude of the magnetic held down¬ 
stream of the shock. For the magnetostatic model, the perpen¬ 
dicular components are scaled manually by the compression ra¬ 
tio, whereas the theory of Schlickeiser (2002) for Alfven wave 
transmission through parallel shocks leads to a more than twice 
as strong downstream magnetic held. As such, with both models 
leading to comparable results for the power spectral indices, one 
might argue that, despite the increased complexity of the Alfven- 
wave model of turbulence is unnecessary. However, the general 
behavior is different with regard to the evolution of the coupled 
wave-particle system. In the presence of self-generated Alfven 
waves, the relative contribution to the overall acceleration of the 
shock and the wave held can vary and lead to a wrong estimation 
of shock parameters if one were to neglect the influence of these 
waves. This effect may therefore complicate the interpretation 
of astrophysical gamma radiation in terms of shock acceleration 
(Aharonian et ak, 2004). This subject is currently under active 
investigation, and results will be presented in due time. In addi¬ 
tion, note that the model used here is strongly restricted by the 
underlying theoretical understanding of Alfven-wave transmis¬ 
sion, which is so far only valid for parallel shock configurations. 

In addition, the work presented here investigated the shock 
acceleration process on a non-relativistic particle ensemble with 
a Maxwell-Boltzmann distributed rigidity injection profile. In 
conformance with the findings of Caprioli & Spitkovsky (2014), 
power spectra with spectral indices of close to -1.5 are obtained, 
and no injection problem is apparent. Similar to the relativistic 
simulations, which used a constant rigidity injection profile, the 
acceleration process becomes less efficient for oblique shocks. 
For a perpendicular shock then, a spectral index of ca. -2 is 
found. Caprioli & Spitkovsky (2014), too, reported a significant 
decrease of the acceleration efficiency for quasi-perpendicular 
shocks with angles above 45°. Their results also show that paral¬ 
lel shocks are most efficient at accelerating particles, compatible 
to our findings. But, for quasi-parallel shocks, they also report a 
drop in acceleration efficiency, which cannot be seen in the re¬ 
sults presented here, where no significant change in the spectral 
index is observable for oblique shocks below 45°. 

A limitation of the results presented here is the slab geometry 
that was used in the magnetostatic model of turbulence applied 
throughout the simulations presented in this work. This was cho¬ 
sen to ensure both the applicability of the Alfven wave transmis¬ 
sion and reflection and the comparability with the magnetostatic 
model. However, a more realistic model might have a poten¬ 
tially large influence on the quasi-perpendicular shock simula¬ 
tions, because this geometry is symmetric in the directions per¬ 
pendicular to the magnetic background field. Thus, in the limit of 
a perpendicular shock, the total turbulent magnetic held is static, 
even when it flows with the ambient plasma upstream or down¬ 
stream of the shock. As such, particles only see the increase of 
the magnetic held strength, but are otherwise unaffected by the 
shock front. Future work certainly has to conduct the simula¬ 
tions presented here for different turbulence geometries, in or¬ 
der to investigate the effect this has on the acceleration process. 
Potentially, a more sophisticated geometry, such as the one pro- 
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posed by Goldreich & Sridhar (1995), will increase the diffusiv- 
ity along the shock normal for quasi-perpendicular shocks and 
thereby improve the acceleration efficiency, increasing the num¬ 
ber of high-energy particles and yielding higher spectral indices 
in the energy spectrum of the particles. 

In conclusion, it has been shown that, even with the lim¬ 
ited configurations used by the test-particle simulation presented 
in this work, reliable results can be obtained. This lays the 
necessary foundation for future studies. Most notably, the ef¬ 
fect of the turbulence geometry on the acceleration efficiency 
of quasi-perpendicular shocks should be investigated, as men¬ 
tioned earlier. Additionally, the influence of many other simu¬ 
lation parameters has to be checked. These parameters include, 
but are not limited to, the compression ratio, the shock speed, 
as well as the energy spectrum of the magnetic turbulence and 
the relative turbulence strength. Furthermore, if an extension of 
the Alfven wave transmission theory of Schlickeiser (2002) to 
oblique shocks can be derived, then it should be tested here 
as well. Also, a more sophisticated model of the background 
magnetic field, e. g., following the research by Uyaniker et al. 
(2002); Eriksen et al. (2011), could be included. Finally, one 
should lessen the simplifying assumptions on the shock front 
itself. Instead, a three-dimensional model of an expanding shock 
wave of finite thickness, potentially including precursor instabil¬ 
ities and/or Rayleigh-Taylor instabilities, could be used. In such 
a model, the shock front will have a varying obliqueness toward 
the magnetic background field, with so-far unexplored effects on 
the final upstream particle energy spectrum. 
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